Field Propagation using Runge-Kutta integration and example15 as demo #166

jonapost · 2021-12-14T18:13:37Z

Development of Runge-Kutta integration in 'kernel' code and example15 to demonstrate its use.

New classes that implement the RK integration include:

RkIntegrationDriver
DormandPrince45
MagneticFieldEquation
fieldPropagatorRK

The development is not yet fully mature.

Runge-Kutta integration has been checked in the unit test vs helix for a uniform field. The classes to integrate it into this example are created, and initial tests done.

phsft-bot · 2021-12-14T18:16:02Z

Can one of the admins verify this patch?

magneticfield/inc/fieldPropagatorConstBz.h

hahnjo · 2021-12-15T10:44:22Z

magneticfield/inc/fieldPropagatorConstBz.h

+        // Also - should we just re-use endPosition / endDirection (in place of midPosition etc ) ??
+
+        // Check
+        accurateIntersection = (midPosition - position).Mag2() < deltaIntersectSq;


Shouldn't this check (midPosition - intersectionPos)? I think at this point position is still the start point of the chord step.

Yes, you are correct! Well spotted. I will fix it.

hahnjo · 2021-12-15T16:08:39Z

I tried to do some testing on my end, here are the results:

A full build of this PR / branch fails, as also reported by Jenkins:

In file included from /home/jhahnfel/AdePT/src-rk-integration/test/test_magfieldRK.cpp:14:
/home/jhahnfel/AdePT/src-rk-integration/magneticfield/inc/PrintFieldVectors.h: In instantiation of ‘void PrintFieldVectors::PrintSixvecAndDyDx(const Real_t*, int, const Real_t*, const Real_t*) [with Real_t = float]’:
/home/jhahnfel/AdePT/src-rk-integration/test/test_magfieldRK.cpp:86:40:   required from ‘bool TestEquation(const Field_t&) [with Real_t = float; Field_t = UniformMagneticField]’
/home/jhahnfel/AdePT/src-rk-integration/test/test_magfieldRK.cpp:457:83:   required from here
/home/jhahnfel/AdePT/src-rk-integration/magneticfield/inc/PrintFieldVectors.h:106:12: error: no matching function for call to ‘Print3vec(const float&, const char [25])’
   Print3vec( dydx[3], " dy/dx [3-5] (=dP/ds) = " );  printf("\n");
   ~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

example15 runs (with quite verbose output), but the simulated results for 10 GeV electrons in CMS with B_z = 1T (all default values) look very wrong:

Mean energy deposit          0.247005 GeV
Mean number of gamma         90.23
Mean number of e-            301.51
Mean number of e+            4.96
Mean number of charged steps 1471.25
Mean number of neutral steps 10911.3
Mean number of hits          11318.4
Killed in propagation:       675

compare that to example13 with also B_z = 1T enabled:

Mean energy deposit          9.87802 GeV
Mean number of gamma         3563.5
Mean number of e-            11807.4
Mean number of e+            309.65
Mean number of charged steps 29275.1
Mean number of neutral steps 32100.3
Mean number of hits          26898.8
Killed in propagation:       392

I tried increasing max_iterations = 1000, but not a huge difference.

The changes to the helix propagator driver to find more accurate intersections yield different results in TestEm3, but still not very good: In particular the first layers have a much lower energy deposit than master and Geant4, going as low as 2 per-mille relative error for the second layer. The later layers have a higher mean energy deposit (compared to master or my attempts to fix it), but still 0.05% below G4HepEm on the CPU. This is probably caused by less energy being deposited in the first few layers and therefore more energetic particles reaching the rest of the detector and "correcting" the underestimated energy deposit.

jonapost · 2021-12-15T21:04:28Z

Thanks for the feedback.
In the last day or two I am having difficulties using cuda-gdb on titanx to investigate. This has slowed me down.
Unless it works I need to move to (or use in parallel) a different machine where I can use it.

hahnjo · 2022-01-31T09:08:42Z

FWIW the very low energy deposit I reported in #166 (comment) is still there.

jonapost · 2022-02-22T11:59:42Z

I have fixed a key bug in the RK glue code for GPU. I get almost perfect agreement between helix and Runge-Kutta now.
There is some cleanup to do, deleting the check code, before running further tests.

hahnjo

I gave this another try, and after fixing the build failure related to kB2C in fieldPropagatorConstBz, it's at least running. However, the energy deposit is still very wrong: With stepLengthFromPhysics = 0.1, I see around 0.5 GeV for a 10 GeV primary. With that line removed, I get only 3-4 MeV (!).

I also think the CompareResponseVector3D isn't suitable for comparing directions, which are / should be normalized so their magnitude is always (close to) 1.

hahnjo · 2022-02-23T08:07:16Z

examples/Example15/electrons.cu

+
+    float BzFieldValue= *gPtrBzFieldValue_dev;   // Use vecgeom::Precision ?
+    if (BzFieldValue != 0.0) {
+      stepLengthFromPhysics= 0.1;


I suppose this is from debugging?

Yes, of course. I am making revisions - I had hoped to finish and push them already yesterday, but got sidelined.

jonapost · 2022-02-24T16:37:26Z

I confirm that I have problems on steps which require multiple integration steps at least.
I am investigating.

I note that the function that compared vectors checked the change from the previous to the next vector - not the magnitude of each vector. So it was not a complete check, but it was a valid check.

jonapost · 2022-02-26T07:09:27Z

I am still investigating the cause of the errors. They occur when multiple steps are needed to complete an integration.

I pushed the code including debugging prints, so that it is in this repository.

hahnjo · 2022-11-10T10:59:21Z

It appears this example is just generally broken, even without magnetic field: Both for trackML.gdml and cms2018.gdml, there is a energy deposit of only 5 MeV, no matter the energy set for the primaries (tested from 1 GeV to 100 GeV). For testEm3.gdml, the energy deposit seems about right (in agreement with example13).

agheata · 2022-11-10T11:05:34Z

It appears this example is just generally broken, even without magnetic field: Both for trackML.gdml and cms2018.gdml, there is a energy deposit of only 5 MeV, no matter the energy set for the primaries (tested from 1 GeV to 100 GeV). For testEm3.gdml, the energy deposit seems about right (in agreement with example13).

But wasn't this just copied from one of the previous examples replacing just the field propagator? @jonapost ?

magneticfield/inc/fieldPropagatorConstBz.h

jonapost · 2022-11-24T17:40:42Z

In example15 I am getting energy deposition which is too perfect when I run with cms2018.gml - a 10 GeV primary deposits exactly 10 GeV, and a 5 GeV primary exactly 5 GeV.

Using 10,000 primaries of 10 GeV the results are:

Mean energy deposit          10 GeV
Mean number of gamma         3601.24
Mean number of e-            11665.6
Mean number of e+            315.006
Mean number of charged steps 24934.3
Mean number of neutral steps 33835.6
Mean number of hits          25198.5

and with 10,000 primaries of 5 GeV:

Mean energy deposit          5 GeV
Mean number of gamma         1805.51
Mean number of e-            5559.16
Mean number of e+            162.994
Mean number of charged steps 13560
Mean number of neutral steps 20526.9
Mean number of hits          17244.2

hahnjo · 2022-11-25T08:31:34Z

Ooohhhh, I think I see my mistake: All other recent examples default to cms2018.gdml, while example15 uses trackML.gdml. That's a bit unfortunate because it is mostly filled with Air and only a bit of Silicon, so there is only very little physics happening. Can we change this back?

jonapost · 2022-11-25T12:07:53Z

Yes, I agree and will change it back. I had changed it for my testing, but it is simpler and better to have the same default.

jonapost · 2022-11-25T16:25:58Z

I am working on the rebase. It is a complicated because of the changes introduced in example13, from the use of the RNG to more local variables in place of data structures ...

agheata · 2023-01-13T09:07:50Z

I am working on the rebase. It is a complicated because of the changes introduced in example13, from the use of the RNG to more local variables in place of data structures ...

@jonapost any progress with the rebase?

agheata · 2023-10-17T07:49:19Z

@jonapost looking at the modified code, I see that it is well-localized in magneticfield/inc and example15, example13 is barely touched. I can help with rebasing this branch, or if it turns out difficult we can make a new one that just adds the relevant files.

agheata · 2023-10-17T08:08:37Z

Actually @jonapost I did the simple exercise and checked out your branch and rebase on the master, which gave no conflict. Then I could even make it compile by changing fieldPropagatorConstBz.h which seems to contain some debugging code compared to the master.

-    vecgeom::NavStateIndex &next_state, bool &propagated, const Precision safety,  int indx, // Slot index 
+    vecgeom::NavStateIndex &next_state, bool &propagated, Precision safety, int indx, // Slot index

- MagneticFieldEquation - DormandPrince45 - Stepper - RkIntegrationDriver - fieldConstants.h with kB2C, delta_intersection parameters - fieldPropagatorRungeKutta : new class for use by electrons.cu, used in new Example14 More: Example15: field propagation using Runge-Kutta. Based on Example13 otherwise. Tests - unit test test_magfieldRK.cpp Simple checks for equation, stepper and driver classes Driver: check version of Advance and V1 (old) Checks driver vs helix results (on cpu). Details: -------- RkIntegrationDriver Introduced new Advance method. Using AdvanceV1 as backward compatible. Enabled used of RK classes with double integrands (field remains float.) PrintFieldVectors: auxiliary methods, initially host-only Changed VECCORE_ATT_HOST_DEVICE to __host__ __device__ example15: Added ability to print Track info check result of RK integration in electrons.cu report differences Optional argument for Bz field value. Use TrackML as default geometry.

…n printing.)

…s RK.

… flip side if failing After a maximum number of failed (zero) steps, reducing the step size each time, accept that the track will not cross the boundary, and flip it to the other volume, which it is apparently entering (or never left.)

…field/inc/CompareResponses

…fety in args

boundary; fieldPropagatorConstBz::ComputeStepAndNextVolume - Optional defaul values in interface (at compile time) - Protected some verbosity using if(verbose) - Added some optional verbosity Both of these are to be trimmed / refined.

…puteNextStepAndVolume Enable defaults for the last 3 arguments of ComputeNextStepAndVolume

…arg is const

Safety was const in ComputeStepAndNextVolume implementation only - added it to declaration.

…rop...

jonapost · 2023-12-04T09:35:22Z

I can run it now successfully with the CMS 2018 geometry for 100,000 initial particles.
So I think it's ready to merge. Can you please clarify what's wrong, if anything?

agheata · 2023-12-04T10:22:34Z

I can run it now successfully with the CMS 2018 geometry for 100,000 initial particles. So I think it's ready to merge. Can you please clarify what's wrong, if anything?

Thanks @jonapost , do you have an estimate timing RK versus uniform field for CMS?

agheata · 2023-12-04T10:26:53Z

The CI shows the example15 failing at the moment: https://lcgapp-services.cern.ch/spi-jenkins/job/AdePT-CI/510/console

jonapost · 2023-12-04T16:35:48Z

The CI shows the example15 failing at the moment: https://lcgapp-services.cern.ch/spi-jenkins/job/AdePT-CI/510/console

Oops, I forgot to commit the change of Stack Limit to 8K. Now done.

jonapost · 2023-12-04T17:10:57Z

My reading is that it was successful : https://lcgapp-services.cern.ch/spi-jenkins/job/AdePT-CI/511/console

jonapost · 2023-12-04T17:21:16Z

I can run it now successfully with the CMS 2018 geometry for 100,000 initial particles. [....]

[..], do you have an estimate timing RK versus uniform field for CMS?

I am adding the numbers is overleaf:

For CMS-2018 and 10,000 electrons (of 10GeV each) the runtime on omenrtx in seconds are:

	B= 1.0 T	B= 3.8 T	Example Name
helix	235.9(3)	571.4(5)	ex13
helix/short steps	231.0	580.5(9)	ex18
Runge-Kutta	413.0	1296 (1)	ex15

agheata · 2023-12-04T17:37:23Z

@Japost, I also measured that Exa15 is slower than Exa13 with a factor of 2 for smaller field values and a bit less for 3.8T. I think we can merge this, but we should be careful in quoting the results from this early version before we understand better the main bottlenecks and if we can easily remove some.

agheata · 2023-11-16T16:50:20Z

examples/Example15/README.md

+SPDX-License-Identifier: CC-BY-4.0
+-->
+
+## Example 14


Suggested change

## Example 14

## Example 15

jonapost changed the title ~~Japost/rk integration~~ [Work in progress] Field Propagation using Runge-Kutta integration and example15 as demo Dec 14, 2021

jonapost changed the title ~~[Work in progress] Field Propagation using Runge-Kutta integration and example15 as demo~~ Field Propagation using Runge-Kutta integration and example15 as demo [Work in progress] Dec 14, 2021

jonapost requested a review from hahnjo December 14, 2021 18:14

jonapost requested a review from agheata December 14, 2021 18:21

hahnjo reviewed Dec 15, 2021

View reviewed changes

jonapost force-pushed the japost/rk-integration branch 2 times, most recently from 61d6bd8 to c6d90b3 Compare February 22, 2022 11:57

hahnjo reviewed Feb 23, 2022

View reviewed changes

jonapost force-pushed the japost/rk-integration branch from fbdb0ba to 3b2c27c Compare June 2, 2022 15:02

hahnjo reviewed Nov 17, 2022

View reviewed changes

magneticfield/inc/fieldPropagatorConstBz.h Outdated Show resolved Hide resolved

WitekPokorski mentioned this pull request Nov 21, 2022

Support for non-constant magnetic field #150

Open

jonapost force-pushed the japost/rk-integration branch from f947388 to 7325e87 Compare February 10, 2023 14:59

jonapost force-pushed the japost/rk-integration branch from 7325e87 to f3ccb0b Compare October 25, 2023 11:39

jonapost added 13 commits November 8, 2023 11:32

example15.cu: fix call to ensure that ClearQueue occurs (not just whe…

3be869b

…n printing.)

example15/electrons.cu: first cleanup, removed max_step=0.25mm, alway…

9713b67

…s RK.

example15.cu: Added optional printing and info printing ; Revised mag…

fc01ae8

…field/inc/CompareResponses

example15: reset cms2018.gdml as default geometryf

fdae177

example15.cu: removed most (optional) verbosity

c36361b

fieldPropagatorConstBz::ComputeNextStepAndVolume: moved slot after sa…

114d404

…fety in args

Example13/electrons.cu : added arguments to ComputeStepAndNextVolume

ee1591f

fieldPropagatorConstBz: cleaned up several printouts; defaults in Com…

1c5a8b5

…puteNextStepAndVolume Enable defaults for the last 3 arguments of ComputeNextStepAndVolume

fieldPropagatorRungeKutta: removed optional prints & cleanup

4787f69

fieldPropagatorConstBz::ComputeStepAndNextVolume enforce that safety …

4225101

…arg is const

fieldPropagatorConstBz.h : clean-up; added const in arg

1b7f554

Safety was const in ComputeStepAndNextVolume implementation only - added it to declaration.

jonapost force-pushed the japost/rk-integration branch from 604fdc7 to 1b7f554 Compare November 8, 2023 10:33

jonapost added 5 commits November 16, 2023 16:34

fieldPropagatorConstBz: fix vecgeom::Precision; cosmetics

d033886

Example13: set B=3.8T ; electrons.cu: took out slot in args of fieldP…

efa21e6

…rop...

CompareResponces.h: Added Copyright, reuse

8c11d26

Example15/electrons.cu : cleanup. Fix to remove slot.

0d2cc3a

example15.cu: cleanup

edbb17a

jonapost changed the title ~~Field Propagation using Runge-Kutta integration and example15 as demo [Work in progress]~~ Field Propagation using Runge-Kutta integration and example15 as demo Dec 4, 2023

example15.cpp: Set stack limit to 8K

e4cccd5

agheata approved these changes Dec 5, 2023

View reviewed changes

agheata merged commit 8e7b9b4 into apt-sim:master Dec 5, 2023
2 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Field Propagation using Runge-Kutta integration and example15 as demo #166

Field Propagation using Runge-Kutta integration and example15 as demo #166

jonapost commented Dec 14, 2021

phsft-bot commented Dec 14, 2021

hahnjo Dec 15, 2021

jonapost Jan 5, 2022

hahnjo commented Dec 15, 2021

jonapost commented Dec 15, 2021

hahnjo commented Jan 31, 2022

jonapost commented Feb 22, 2022

hahnjo left a comment

hahnjo Feb 23, 2022

jonapost Feb 24, 2022

jonapost commented Feb 24, 2022

jonapost commented Feb 26, 2022

hahnjo commented Nov 10, 2022

agheata commented Nov 10, 2022

jonapost commented Nov 24, 2022 •

edited

Loading

hahnjo commented Nov 25, 2022

jonapost commented Nov 25, 2022

jonapost commented Nov 25, 2022

agheata commented Jan 13, 2023

agheata commented Oct 17, 2023

agheata commented Oct 17, 2023

jonapost commented Dec 4, 2023

agheata commented Dec 4, 2023

agheata commented Dec 4, 2023

jonapost commented Dec 4, 2023

jonapost commented Dec 4, 2023

jonapost commented Dec 4, 2023 •

edited

Loading

agheata commented Dec 4, 2023

agheata Nov 16, 2023

Field Propagation using Runge-Kutta integration and example15 as demo #166

Field Propagation using Runge-Kutta integration and example15 as demo #166

Conversation

jonapost commented Dec 14, 2021

phsft-bot commented Dec 14, 2021

hahnjo Dec 15, 2021

Choose a reason for hiding this comment

jonapost Jan 5, 2022

Choose a reason for hiding this comment

hahnjo commented Dec 15, 2021

jonapost commented Dec 15, 2021

hahnjo commented Jan 31, 2022

jonapost commented Feb 22, 2022

hahnjo left a comment

Choose a reason for hiding this comment

hahnjo Feb 23, 2022

Choose a reason for hiding this comment

jonapost Feb 24, 2022

Choose a reason for hiding this comment

jonapost commented Feb 24, 2022

jonapost commented Feb 26, 2022

hahnjo commented Nov 10, 2022

agheata commented Nov 10, 2022

jonapost commented Nov 24, 2022 • edited Loading

hahnjo commented Nov 25, 2022

jonapost commented Nov 25, 2022

jonapost commented Nov 25, 2022

agheata commented Jan 13, 2023

agheata commented Oct 17, 2023

agheata commented Oct 17, 2023

jonapost commented Dec 4, 2023

agheata commented Dec 4, 2023

agheata commented Dec 4, 2023

jonapost commented Dec 4, 2023

jonapost commented Dec 4, 2023

jonapost commented Dec 4, 2023 • edited Loading

agheata commented Dec 4, 2023

agheata Nov 16, 2023

Choose a reason for hiding this comment

jonapost commented Nov 24, 2022 •

edited

Loading

jonapost commented Dec 4, 2023 •

edited

Loading